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The Putative Liquid-Liquid Transition is a Liquid-Solid Transition in 
Atomistic Models of Water 

David T. Limmer^ and David Chandler-'^'r^l 

Department of Chemistry, University of California, Berkeley, California 94720 

(Dated: 21 September 2011) 

We use numerical simulation to examine the possibility of a reversible liquid-liquid transition in supercooled 
water and related systems. In particular, for two atomistic models of water, we have computed free energies 
as functions of multiple order parameters, where one is density and another distinguishes crystal from liquid. 
For a range of temperatures and pressures, separate free energy basins for liquid and crystal are found, 
conditions of phase coexistence between these phases are demonstrated, and time scales for equilibration are 
determined. We find that at no range of temperatures and pressures is there more than a single liquid basin, 
even at conditions where amorphous behavior is unstable with respect to the crystal. We find a similar result 
for a related model of silicon. This result excludes the possibility of the proposed liquid-liquid critical point 
for the models we have studied. Further, we argue that behaviors others have attributed to a liquid-liquid 
transition in water and related systems are in fact reflections of transitions between liquid and crystal. 

Keywords: water, liquid-liquid transformations, critical behavior, free energy 



I. INTRODUCTION 

This paper reports the results of a numerical study 
aimed at elucidating the purportecP^l liquid- liquid phase 
transition in supercooled liquid water. The results in- 
dicate that this hypothesized polyamorphism does not 
exist in atomistic models of water. While not contradict- 
ing the existence of irreversible polyamorphism of the sort 
observed in non-equilibrium disordered solids of water ,'2n2l 
and not excluding the possibilities of liquid-liquid transi- 
tions in liquid mixturegSl, polymerizing fluidgSl and some 
theoretical models,'I2Hl3 ^j^g results do suggest that a re- 
versible transition and its putative second critical point 
are untenable for one-component liquids, like water, that 
exhibit local tetrahedral order and freeze into crystals 
with similar but extended order. 

The terminology "transition" is used here to refer to 
distinct phases, where coexistence implies the formation 
of interfaces that would spatially separate the coexisting 
phases or to response functions that diverge in the ther- 
modynamic limit .1^ The structural changes for a transi- 
tion between two liquids or between a liquid and a crystal 
are distinct from continuous pressure induced changes in 
normal liquid water. ^^ These changes associated with a 
phase transition are global and therefore are also distinct 
from bi-continuous behaviors that do not persist beyond 
small length scales. ^^ 

Polyamorphism of water has been achieved through 
various out-of-equilibrium experimental protocols result- 
ing in a multitude of thermodynamically unstable, kinet- 
ically trapped structures.'^2"22l These different disordered 
structures have been generally partitioned into two gen- 
eral categories known as either low-density amorphous 
gQiy^UHlSl or high-density amorphous solids. '22H221 Some 



'Electronic mail: cliandler@berkeley.edu 



have interpreted changes in these structural motifs as 
non-equilibrium manifestations of an underlying equilib- 
rium phase transition between two forms of liquid water. 
This conjecture forms the basis of some attempts to ex- 
plain many of the well-known anomalous thermodynamic 
properties of water, e.g., Refs. [T| [23l [24l and [25l There 
are other ways of explaining these anomalies ,123 but the 
phase transition hypothesis seems particularly intriguing, 
and it is the focus of this paper. 

The hypothesis is impossible to test by natural experi- 
ments because the location of the presumed transition is 
outside experimentally accessible conditions.'^ In partic- 
ular, bulk supercooled water is unstable as a liquid, and it 
rapidly crystallizes in the regime of predicted polyamor- 
phism. While the properties of non-equilibrium glassy 
materials can be studied in this region, it is uncertain 
whether inferences regarding reversible thermodynamic 
behavior can be made from such measurements. Some 
experiments have studied water confined to long pores 
with radii no larger than 1 nm.'^SHSl] These experiments 
avoid the instability and thereby attempt to detect man- 
ifestations of the transition. While water does not freeze 
in such circumstances, it is questionable whether proper- 
ties of bulk water can be inferr ed from behaviors of these 
one-dimensional systems.l32E3 

Molecular simulation provides a means to overcome 
this ambiguity. Specifically, sufficiently realistic models 
can be studied computationally while controlling order 
parameters that distinguish liquid from crystal. It is 
in this way that we examine the reversible behavior of 
models of water. Along with establishing coexistence be- 
tween liquid and crystal, we are able to study the dy- 
namics of the transition between these phases. We also 
locate and explore the free energy surface for the region 
of the pressure-temperature phase diagram known as "no 
man's land".^ This is the region where amorphous be- 
havior would be unstable in the absence of control. Our 
results indicate that some observations attributed by oth- 



ers as manifestations of a liquid-liquid transition are in 
fact observations of the temperature-pressure boundary 
separating the region of amorphous instability from that 
of a single phase of amorphous metastability. 

Others have used molecular simulation for realistic 
models of water^^^ and related liquida^^SSl iq exam- 
ine their possible polyamorphisms. Those cited here^BfSl 
are representative but by no means comprehensive (we 
exclude from this list models that do not exhibit local 
tetrahedral order, e.g. Ref. 10). In all cases, the meth- 
ods employed have been limited in at least one of three 
ways: time scales that are too short, system sizes that are 
too small, and order parameters that fail to discriminate 
order from disorder or fail to be adequately controlled. 
By employing multiple order parameters and free-energy 
sampling methods, we overcome time-scale issues and are 
able to discriminate between phases of different symme- 
tries. By considering different system sizes and size scal- 
ing analysis, we overcome uncertainty associated with fi- 
nite system sizes. 

Most of the results we present in this paper have been 
computed with a recently developed model by Molinero, 
so-called "mW" water .HS We use this model for three rea- 
sons. First, it is a computationally convenient model be- 
cause it contains no long-ranged forces, relying instead 
on short-ranged three-body forces to favor microscopic 
structures consistent with those of water. Second, the 
behavior of the model is realistic in the sense that in the 
range of conditions we wish to study its phase diagram 
is a reasonable caricature of that for water.^'*^ ^^ Third, 
the results obtained with this model would seem to apply 
to other systems in addition to water in that the model 
is a variant of one developed by Stillinger a nd Webe r,B3 
which has been used to treat behaviors of sPSHlIMl and 
SiOa.SSl 



II. PHASE SPACE STUDIED 

The mW modeP^J differs from the Stillinger- Weber 
inter-particle potential energy function for silicon in two 
ways. The first is the adoption of different values for the 
length and energy parameters of the model. This differ- 
ence is inconsequential to our study because it amounts 
to a simple rescaling of temperature and density. The 
second is more substantive but slight. Specifically, to 
capture some thermodynamic properties of water, the 
mW model has a partitioning between two- and three- 
body terms that differs by 10% from the partitioning for 
silicon. See Ref. |44]for details concerning the niW model. 



A. Pressure-temperature phase diagram of the mW model 

We have studied the phase behavior of the mW model 
by computing free energy surfaces throughout its con- 
densed phases. Figure IT] shows the state points exam- 
ined. Each circle represents a state point where the free 
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FIG. 1. Phase space sampled in our calculations for the mW 
model. Open circles refer to states where the amorphous liq- 
uid is found to be unstable, filled red circles refer to states 
where the amorphous liquid is found to be metastable with 
respect to the crystal, and filled blue circles refer to states 
where the liquid is found to be stable with respect to the crys- 
tal. The labeled circles (a), (b) and (c) identify state points 
where expHcit free energy surfaces are shown in Fig. [S] The 
circles labeled (d) and (e) identify state points where explicit 
free energy surfaces are shown in Fig. HI The black circle 
labeled (f) is the state point where crystal-liquid phase coex- 
istence is examined in Fig. [5] Lines, diamonds and squares 
locate previous estimates of a liqu id-liquid phase transition 
inferred from experimental results p^'^^'^ see text. A star 
locates a previous prediction of a liquid-liquid critical point 
based upon extrapolation of simulation results for the mW 
model.HSl The blue triangle and hexagon are estimates of low 
temperature critical point locations obtained from interpret- 
ing simulation results for the mSTz^ model and the ST2i^ 
model, respectively. 



energy has been calculated as a function of the global 
system density and an order parameter that quantifies 
broken orientational symmetry. The latter distinguishes 
an amorphous liquid phase from a crystal, whereas the 
former distinguishes two amorphous phases with differ- 
ent densities. Basins in the free energy surface establish 
relative stabilities of the phases. 

To put this diagram in context, we highlight state 
points that others have identified as relevant to a liquid- 
liquid phase transition in supercooled water. The tem- 



perature of maximum density at low pressure sets the 
scale of the figure. This chosen reference temperature is 
To = 250K for the mW modelP and it is Tq = 277K 
for water !^ The phase diagram in Fig. [l] shows that 
the density maximum of liquid mW occurs at slightly 
supercooled conditions while that of experimental water 
occurs at a temperatures slightly higher than the freezing 
temperature. 

The points identified by Liu et al.^* come from 
measured relaxation times of water confined in silica 
nanopores with a 7A radius. These relaxation times 
have a temperature dependence that changes from super- 
Arrhenius to Arrhenius upon cooling below a crossover 
temperature, T^{p). This temperature depends upon 
external pressure p, and points on this line are shown 
in Fig. [l] with unfilled squares, which are attributed in 
Ref. 28 to crossing a "Widom line" . A Widom line refers 
to a locus of maximum response that ends at a critical 
point.'SZl In cases where a phase transition exists, there 
are many such lines because different response functions 
have different lines of extrema. Ambiguity ceases only in 
the proximity of a critical point. But whether any such 
lines can be related to T^{p) is unclear because Widom 
lines refer to time-independent thermodynamic behavior 
and Tx{p) refers to time-dependent non-equilibrium be- 
havior. 

A different basis for identifying relevant points^ is made 
by Zhang et al. considering the same system.'^''^ In this 
case it is the density of the water that is measured. 
This observed density exhibits hysteresis upon alternat- 
ing heating and cooling scans, and the hysteresis grows 
upon increasing pressure. We have already noted that it 
is questionable whether the phase behavior of bulk water 
can be related to that of water confined to narrow pores. 
Virtually all molecules in those pores are influenced by 
interfaces. Nevertheless, points of maximum hysteresis, 
denoted by the filled diamonds in Fig. [l] have been at- 
tributed to a line of first-order liquid-liquid transitions, 
and points of less significant hysteresis, marked by open 
diamonds in the figure, are attributed to a continuation 
of that transition.ES 

Another proposed line of liquid-liquid transitions is 
constructed by Fuentevilla and Anisimov.l^ Here, a pos- 
tulated scaling form is used to extrapolate from exper- 
imentally accessible equilibrium thermodynamic data. 
The resulting prediction and its analytic continuation 
are drawn as solid and dashed lines in Fig. [l] Even if 
a critical point is present this predicted line is question- 
able because it is generally impossible to identify critical 
divergences from a small rise in noncritical background 
fluctuations of the sort contributing to the heat capacity 
at standard conditions. This fact is illustrated by Moore 
and Molinero's predicted critical point for mW water.SS 
Its location, the star in Fig. [T] is found by extrapolation 
from a small rise in a response function computed at dis- 
tant thermodynamic conditions. We find no evidence for 
a liquid-liquid transition anywhere near this predicted 
critical point. Rather, it and all other estimates pertain- 



ing to a purported liquid-liquid transition lie close to a 
spinodal associated with crystalization. This finding is 
not inconsistent with Moore and Molinero's more recent 
report that the amorphous phase of the mW model seems 
to be forever changing and impossible to equilibrate at a 
point in the phase diagram where liquid-liquid transitions 
have been suggested.!^ 

Two additional marked points in Fig. [T] the blue tri- 
angle and hexagon, refer to other estimates of a location 
for a liquid-liquid critical point. These are estimates ob- 
tained from extrapolating simulation results for variants 
of the ST2 water modeljISS about which we have more to 
say later. 



B. Anomalous thermodynamics of the mW model 

Water exhibits anomalous thermodynamic properties 
at low temperatures, properties that are are non-singular 
but nonetheless unusual. Because these behaviors have 
been proposed as indicators of a liquid-liquid transition, 
it is important to show that the mW model exhibits such 
behaviors. Specifically, we focus on the density maximum 
as a function of temperature, and the relatively large rate 
of increases upon lowering temperature of both isother- 
mal compressibility and isobaric heat capacity.*^ 

Thus, we have used a constant pressure ensemble to 
compute p = N/{V}, kt = {{SVf)/kBT{V} and Cp = 
{{5Hf)/kBT^ for the mW model. Here, iV, V and H de- 
note number of molecules, volume and enthalpy, respec- 
tively; 5V and 5H denote deviations from mean values 
of V and H, respectively; the pointed brackets denote an 
ensemble average; fee is Boltzmann's constant. 

Figure [2] compares our computed results at ambient 
pressure with those found from experimental observation 
of water.-° As in Fig. [l] we use the low pressure point of 
density maximum as our reference state for these com- 
parisons. Figure [2] shows that the qualitative trends and 
magnitude of anomalies of mW water agree with those 
of experimental water. The low-temperature end of the 
displayed graphs occur at the point where the liquid be- 
comes unstable. Down to that temperature, the growths 
of Kt and Cp are notable but modest in size and far from 
the sort of divergent behavior one would ordinarily as- 
sociate with a critical point or phase boundary. At all 
stable and meta-stable liquid phase states we have stud- 
ied, see Fig. [Tl we find similar nonsingular behavior. 

The lower panels of Fig. [2] show that the trends ob- 
served at 1 bar persist to higher pressures in the mW 
model, with the density maximum temperature decreas- 
ing slightly as pressure increases. These trends are con- 
sistent with experiment.'^ 



III. ORDER PARAMETERS AND FREE ENERGIES 

In this section, we define order parameters and present 
free energy functions of those order parameters. 
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FIG. 2. Average thermodynamic properties as a function of temperature. Upper panels compare mW model results with those 
of experiment at p = 1 bar. The filled blue circles are the calculated results for the mW model, where error estimates are 
one standard deviation. The empty black circles are experimental results taken from Ref. '50' Lower panels show pressure 
dependence of the mW model results. 



A. Measures of crystalline order 

We use two types of order parameters. One is bulk 
density, the other quantifies orientational order. For the 
latter, we use Steinhardt, Nelson and Ronchetti's Qq and 
i/)6- For a finite system analyzed with computer simula- 
tion, these variables prove more convenient than Fourier 
the components of the density. They also prove more 
useful than dynamic measures, which cannot distinguish 
liquid from crystal at supercooled conditions, where dif- 
fusion is slow in the liquid due to glassy dynamics and 
nonzero in the crystal due to defect motion. 

Both Qg and ijjq are functions of a projection of the 
density field into averaged spherical harmonic compo- 
nents. To evaluate Qi, for each water molecule j, we 
calculate the set of quantities 
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where the sum is over those nearest 4 neighbors, n^. 
Yl" {(pij , 9ij) is the £,m spherical harmonic function as- 
sociated with of the angular coordinates of the vector 
T^i ~ ^j joining molecules i and j, measured with respect 



to an arbitrary external frame. Since q^ ^ is defined in 
terms of spherical harmonics, it transforms simply under 
rotations of the system or the arbitrary external frame. 
These quantities are then summed over all particles to 
obtain a global metric 
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and then contracted along the m axis to produce a pa- 
rameter that is invariant with respect to the orientation 
of the arbitrary external frame. 
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The other orientation order parameter we consider, tpi, 
is evaluated by first defining bond variables through lo- 
cal contractions of the qt^my which are reference frame 
independent. 
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and then summing over all of the bonds made between 
molecule i and its nearest 4 neighbors, 
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(5) 
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Finally, the global parameter is obtained by summing 
over all molecules. 
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The mean or most probable value of Qe for an amor- 
phous phase approaches zero in the thermodynamic limit, 
while it is finite for a crystalline phase. As such, Qi 
is a distinguishing order parameter for amorphous and 
crystalline phases. In contrast, because its contractions 
occur locally and not over the entire system, ipi is non- 
vanishing in the thermodynamic limit for both disordered 
and ordered states. Nevertheless it is a useful measure 
of orientational order because the distributions of V'^ for 
the low temperature liquid are sensitive to the amount 
of crystallization in the system, and their mean values at 
low temperatures differ significantly between liquid and 
crystal. Further, as t/ji retains local information, it is use- 
ful in determining the existence of grain boundaries and 
defects. We have taken the £ = 6 multipole because we 
have found empirically that it is particularly sensitive to 
distinguishing liquid water and ice. 



B. Free energy surfaces at conditions of metastability 

Free energies of density p, orientational order parame- 
ters Qe and ipe, and so forth, are related to the probabili- 
ties of the order parameters in the usual way. Specifically, 



F{p,Qe,...)^-kBTliiP{p,Qe 



const. 



(7) 



where the probability P {p,Qe, ■■■) is proportional to the 
partition function for micro-states with the specified val- 
ues of the order parameters. The irrelevant additive con- 
stant in Eq. [7] refers to normalization and standard state 
conventions. 

To evaluate the probabilities and their associated free 
energies, we have adopted a hybrid Monte Carlo simu- 
lation approach as used by Duane et. al.l^We consider 
ensembles with N , p, and T fixed. Two different moves 
are made within this framework: random changes in vol- 
ume, and short molecular dynamics trajectories. These 
are made with a ratio of 1:5. Maximum volume dis- 
placement and maximum molecular dynamics trajectory 
length are adjusted to yield a 30% acceptance. This tech- 
nique produces suitably swift equilibration even within 
the supercooled regime. 

The order parameters p, Qe or ipe, are controlled with 
umbrella sampling, by propagating the system under 
its unbiased hamiltonian and computing the order pa- 
rameters only when determining Metropolis acceptance 



probabilities. All molecular dynamics propagation was 
done using the LAMMPS molecular dynamics simula- 
tion package.!^ Most of the free energy calculations were 
accomplished with 216 particles. For each window in the 
umbrella sampling, the simulations ran long enough to 
obtain at least 1000 independent samples of each of the 
biased observables. The umbrella biasing potentials em- 
ployed were 



AU = k{p-p*)^ + KiQe-Q*ef 



(8) 



or the same formula with Qe replaced by ipQ. Adopting 
K in the range of 500 to 2000 fceT and k in the range 
of 1000 to 2000 ksT cm^/g proved satisfactory. Statis- 
tics gathered in these biased ensembles were unweighted 
and the free energy differences between each ensemble 
were estimated using multi-state Bennett acceptance ra- 
tio (MBAR).^^ Error estimates for the free energies we 
have calculated in this way are less than k^T. 

Figure [3] depicts representative free energies for three 
different state points. The locations of those state points 
are noted in Fig. [l] Each free energy surface includes 
the range of densities where liquid and crystal basins are 
located. With the variable Qe, we see a significant sep- 
aration between liquid and crystal basins. For the state 
points considered in Fig. |3J with N = 216, the crystal 
basin is centered around Qe ~ 0.5, while the liquid basin, 
when it exists, is centered around Qe ss 0.05. As N in- 
creases, the former changes little, but the latter tends to 
zero. This behavior is illustrated explicitly in the next 
section. 

The state points considered in Fig. [3] show how the 
free energy surfaces evolve as the pressure or temperature 
are changed. In Fig. |3]^a) , a barrier separates the liquid 
phase from the crystal. Therefore, at that state point 
the liquid is metastable. Lowering the temperature and 
increasing the pressure. Fig. p^ (b) shows the barrier to 
crystallization has vanished. Further decreasing temper- 
ature. Fig. ^ (c) shows increasing driving force towards 
the stable crystal. At those state points, the liquid is 
unstable. 

Similar behavior is found with the free energy of p and 
ipe- This function, F{p,^e): is shown in Fig. H at two 
different state points. Figure Ilia) shows this free energy 
at a temperature and pressure, labeled (d) in Fig. nj 
where the liquid is metastable with respect to the crys- 
tal. At this state point, the mean value, ("06) is about 
0.27, a value that reflects the relatively small amount of 
local ordering present in the supercooled liquid. In con- 
trast, for the crystal we find {tpe) ~ 0.9. Fig. l4Fb) shows 
the free energy for a temperature and pressure in the re- 
gion of the phase diagram where the amorphous phase 
is unstable, the so-called no man's land. This point, la- 
beled (e) in Fig. Ill is close to a proposed location of a 
liquid-liquid critical point.'^ We see, however, that it is 
not a point of criticality. The behavior of ipe is strongly 
correlated to the potential energy. This fact follows from 
the functional form of tpe a-nd the three-body potential 
of the mW model. Thus, the behavior of F{p, ipe) should 
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FIG. 3. Free energy surfaces for mW water as a function of p and Qq. The system is periodically replicated and contains 
A'' = 216 particles. As shown in (a), the liquid is metastable with respect to the crystal. As the system is cooled, the barrier 
disappears, as illustrated in (b). Finally in (c) the free energy obtains a large gradient along the Qq direction and fluctuations 
in density are damped out. Adjacent contour lines are spaced by 1 A:bT, and statistical uncertainties are smaller than that 
energy. 



be similar to that of F{p^ U) where U denotes the total 
potential energy of the mW model. 

For all of the state points considered, which includes 
a broad swath of no-man's land, there is no evidence of 
a bifurcation of the free energy along the density direc- 
tion within the liquid region (i.e., where Qg and V'e s-^e 
small). What bifurcation does exist is associated with 
a transition between an amorphous phase and a crystal. 
We now consider whether it is a first order transition. 
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FIG. 4. Free energy for mW wateras a function of p vs ipa 
at conditions where the liquid is (a) metastable, and (b) un- 
stable. The system is periodically replicated and contains 
A'^ = 216 particles. Adjacent contour lines are spaced by 1 
k^T, and statistical uncertainties are smaller than that en- 
ergy. 



IV. FREEZING TRANSITION 



A. mW Model 



The character of the transition can be analyzed by 
studying the system-size dependence of the contracted 
free energy.E3 



F{Qq) = -fceTln (j dp exp[-/?F(p,Qe 



(9) 







This function is shown in Fig. [5]for the mW model at one 
of the pressures and temperatures where an amorphous 
phase is in coexistence with the crystal. Such points are 
at the boundary between the blue and red regions in Fig. 
T] The quantity AF{Qe) = F{Qe) - min[F{Qe)] reaches 
its maximum value when an interface separating amor- 
phous and crystal phases extends across the entire sys- 
tem. This maximum value is the interfacial free energy. 
Accordingly, for a first-order transition, it should be pro- 
portional to N"^/^. This scaling is satisfied to a good 
approximation for the system sizes considered in Fig. [5] 

Nonzero values of Qq in an amorphous phase are due to 
fluctuations. As such, the mean value of Qg for the amor- 
phous phase should disappear as X/N^/"^ . This scaling is 
also found for the system sizes studied and is illustrated 
in Fig. [5] In contrast, for a crystal Qq will have a nonzero 
mean that remains finite as A'^ — >■ cxi. This behavior is 
consistent with our numerical results, as also illustrated 
in Fig. [5] 

Thus, the transition between liquid and crystal in mW 
water appears to be a standard freezing transition which 
is first order and between phases with different orienta- 
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FIG. 5. Free energies as a function of Qg for A*' — 216 (blue), 
512 (green), and 1000 (red), calculated at r/ro=1.09 and 
p —1 bar. Left inset: Interfacial free energy for different sys- 
tem sizes. For comparison a line of slope 2/3 is also shown. 
Right inset: The mean value of Qe for liquid (circles) and 
crystal (squares) for different system sizes. Error estimates 
are shown in the main figure, but are smaller than the sym- 
bols in the insets. 



tional symmetry. The analysis used here to reach that 
conclusion can be applied to other models. For example, 
we have carried out this analysis to study the Stillinger- 
Weber model of silicon. Here too, we find that the 
model exhibits a freezing transition, and contrary to re- 
cent suggestions there is no evidence for an equilibrium 
liquid-liquid transition. We also arrive at this same con- 
clusion for another model of water, which we turn to now. 



B. mST2 model 

We have considered molecules interacting by a mod- 
ified form of Stillinger and Rahman's pair potential.'^ 
The modification incorporates long-ranged electrostatics 
rather than the simple spherical truncation of the original 
model. The resulting system, which we call "mST2 wa- 
ter," has been studied by Liu et al.,^who report finding 
evidence of a liquid-liquid transition and an associated 
critical point. The mST2 model is more difficult to sim- 
ulate than the mW model because the former contains 
long-ranged interactions and the latter does not. As such, 
our investigation of its behavior is more limited than 
those we have preformed for mW water. Nevertheless, 
our investigation seems sufficient to challenge the finding 
of a liquid-liquid transition at the conditions examined 
by Liu et al. It also seems sufficient to discount an as- 
sortment of less direct simulation studies that also r eport 
evidence of a liquid-liquid transition in ST2 water .ISSESEIl 
Indeed, we have considered this particular model of water 
because it is so often examined in publications supporting 
the hypothesized liquid-liquid transition, most recently in 



a paper'^' motivated by preliminary reports of our work. 

Figure [6] shows free energies we have computed for 
the mST2 model using N — 216 molecules. The pro- 
cedures we employed are identical to those used for the 
mW model, except for the technical detail that we mod- 
ify LAMMPS to handle the specific mST2 potential. We 
focus on the region of the p-T plane where Liu et al. 
report bifurcation in the free energy as a function of den- 
sity. In that region, we too find a bifurcation, but not 
between two amorphous phases. The grand canonical 
Monte Carlo simulation method of Ref. |35] is sufficient 
to detect a phase boundary for the liquid, but it can- 
not distinguish liquid from crystal because it does not 
control distinguishing order parameters. In our calcula- 
tions, where both p and Qe are controlled, we find that 
a boundary does in fact exist between liquid and crystal. 
But at the thermodynamic conditions considered by Liu 
et al., there is no evidence of a second liquid basin in the 
free energy F{p,Qe). 

The specific free energies shown are found by first com- 
puting F{p, Qe) from our simulations at T = 235 K and 
p = 2.2 kbar, i.e. we compute F{p,Q(,) = F{p,Q(i;p,T). 
The free energy shows that for this point in the phase 
diagram, the crystal is stable with respect to the liquid. 
A specific state point considered by Liu et al. is at the 
same temperature but a different pressure or chemical 
potential for which the free energy can be reached by a 
shift in chemical potential 



F{p, Qe;T, Ap) = F{p, Qg) - pVAp, 



(10) 



with A/i « 0.55 ksT. Here, Ap is the chemical potential 
relative to that at {T,p)={235 K, 2.2 kbar), and V is the 
average volume of the system at that temperature and 
pressure. A lower value of Ap brings the system to a 
point of coexistence between the liquid and the crystal. 
These free energy surfaces are shown in Panels (b) and 
(c) of Fig. |6] The free energy computed by Liu et al. is 
the contraction 



F{p) = -fcBTln / dQe exp [-(3F{p,Qe) - ^pVAp] 



(11) 
This function is shown in Panel (a) of Fig. [6] 

Like Liu et al, we find a bistable free energy at this 
temperature and for this size system. The locations for 
the minima wc find for F{p) are in good accord with those 
found by Liu et al. But our free energy has a large barrier 
between the two basins, reflecting a finite crystal-liquid 
surface tension, while that reported by Liu et al exhibits 
a small barrier. Liu et al. suggest that their result is 
indicative of a liquid-liquid transition and the proxim- 
ity of a critical point. However, our free energy surface 
shows no such phase transition behavior. There is only a 
crystal-liquid first-order transition. We suggest that the 
Liu et al. result is a non-equilibrium phenomenon, where 
a long molecular dynamics run at constant T —p and ini- 
tiated from their low-density amorphous phase will even- 
tually equilibrate in either the low density crystal or in 
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FIG. 6. Free energies for the mST2 model of water. The system is periodicaUy replicated and contains A'^ = 216 molecules. 
Panel (a) is the contracted F(p). Panels (b) and (c) are the surfaces F{p,Qs). Phase coexistence between amorphous and 
crystal phases occurs at A/i = 0.27kBT, where A/i is the chemical potential relative to that of phase space point (T, p) = (235 
K,2.2 kbar). Adjacent contour lines in (b) and (c) are spaced by 1 ksT and statistical uncertainties are of the order of, or less 
than, that energy. Error bars in (a) are one standard deviation. 



the higher density metastable liquid. The time scale for 
this equilibration is long, as we discuss in the next sec- 
tion. 

Whatever the cause for the Liu et al results, the bi- 
stability cannot be attributed to a liquid-liquid transi- 
tion without also showing that the barrier separating pre- 
sumed liquid-phase basins satisfies the requisite growth 
with N, scaling as N"^/^. This demonstration has not 
been done, and from our results, it is unlikely that it can 
be done. 

Another variant of the ST2 model, considered by Poole 
et al., louses a reaction field approximation to estimate 
the effects of long-ranged forces. We call it the"ST2r" 
model. One expects similar phase behaviors from the 
mST2 and ST2r models."^ Based on an extrapolation 
from the equation of state computed for ST2r model, 
Poole et al. predict the presence of a liquid-liquid tran- 
sition, and the critical point location obtained from that 
estimate is shown in Fig. IT] The density-maximum refer- 
ence temperature for both mST2 and ST2r is Tq « 330 K. 
Poole et al.'s estimate the critical temperature to be 
Tc = 245 K. Our calculations for mST2, shown in Fig. ^ 
are at the lower temperature, T = 235 K. Accordingly, at 
some pressure, we should find bistable liquid behavior if 
indeed a critical point existed at the higher temperature. 
But we find that upon adding Ap V to our computed 
F(p,(36; 2.2 kbar, 235 K), where Ap ^ p - 2.2kbar, no 
second liquid basin can be discerned for any reasonable 
value of p. Therefore, and similar the to behavior found 
with the mW model, extrapolation from the behavior of 
a one-phase system as done in Ref. [551 proves to be a poor 
indicator of a phase transition. 



V. DYNAMIC METASTABILITY 

To arrive at the results of the previous sections, equi- 
libration is achieved with umbrella sampling. Vari- 
ous other reweighting Monte Carlo procedures could be 
ygg(j 157161] Some researchers, however, attempt to learn 
about a possible reversible phase transition in super- 
cooled water through straightforward molecular dynam- 
ics simulation. This approach is limited to cases where 
relaxation is swift compared to computationally feasible 
trajectory lengths, but relaxation associated with phase 
transitions is generally not swift, especially at super- 
cooled conditions. To judge the feasibility of such an 
approach, it is therefore useful to estimate pertinent re- 
laxation times. For supercooled water, there are two im- 
portant classes: times required to nucleate and grow a 
crystal, and times required to reorganize atomic arrange- 
ments in the liquid. We have estimated both with molec- 
ular dynamics of mW water. We use the equilibrium sam- 
pling described in prior sections to prepare initial config- 
urations from which we carry out Newtonian trajecto- 
ries to compute dynamical properties. These trajectories 
evolve with a Nose-Hoover integratoi^^ with a thermo- 
stat time constant of 5 ps and a barostat time constant 
of 5ps. 

At conditions of liquid metastability, where a free en- 
ergy barrier separates liquid and crystal basins, nucle- 
ation is the rate-determining step to form the equilibrium 
phase. For those conditions, we have computed this rate 
constant following a standard Bennett-Chandler proce- 
dure for rare-event sampling.lSIl Specifically, we take Qe 
as the reaction coordinate, so that the rate constant for 
nucleation is fcnuc = f^exp[— F(Qg)/fcBr], where Qq is 
the point of maximum F(Qq) between liquid and crys- 



tal basins, and the prefactor, v, includes the transmission 
coefficient. This prefactor is determined by sampling tra- 
jectories initialized at the top of the free energy barrier, 
i.e., initialized at configurations with Qg = Qg-'^ Other 
choices of transition state are possible, but the net result 
is invariant to that choice.'^ The mean time to nucleate 
the crystal is then l/fc^uc = ''xti- Results obtained in 
that way with N = 216 mW particles are shown in Fig. 
[7j The reference time used to represent these results, tq, 
is the structural relaxation time at the reference liquid 
state used throughout this paper, Tp — 2bQK. For mW 
this time is tq ~ 0.5 ps; for experimental water, this time 
is larger by a factor of about 5. 

For conditions of liquid instability, (i.e. the no-man's 
land where there is no barrier between liquid and crystal) , 
the method of rare-event sampling is no longer appro- 
priate. For those conditions, we compute first-passage 
times.!^ The results obtained depend upon the initial 
preparation of the system because the unstable system 
is far from equilibrium. In the particular preparation we 
employ, we equilibrate the system in the liquid region at 
T/Tq = 0.84 where the liquid is metastable. Then at 
time i = 0, the system is quenched to the target tem- 
perature and allowed to evolve towards the crystal state. 
The first-passage time is taken as the first time a trajec- 
tory with initial conditions prepared in that way reaches 
a configuration with Qq = 0.2. We find an exponen- 
tial distribution of first-passage times. Mean values from 
that distribution are the values of Txti shown in Fig. [7] 
for no-man's land state points. 

The line in the p-T plane separating filled and unfilled 
circles in Fig. 1 is the boundary between metastable 
and unstable liquid conditions. At metastable conditions 
close to that boundary, we have checked that the Txti 
found from the first-passage method agrees with that 
found from the rare-event sampling. 

For the structural relaxation time of the supercooled 
metastable liquid, we consider trajectories initiated from 
equilibrated configurations in the liquid region and ob- 
serve the time t it takes an average particle to move one 
molecular diameter, i.e., 3A=(1/A^) ^^ |ri(i) — ri(0)|. For 
temperatures below the limit of liquid stability, we use 
the same procedure for generating initial conditions as 
we used in the calculation of the first passage times. At 
the higher temperatures we considered, we find an ex- 
ponential distribution of these times. At temperatures 
below T/Tq w 0.88 the distribution deviates from an 
exponential, and increasingly so as temperature is fur- 
ther lowered. This behavior implies the onset of glassy 
dynamics^ with an onset temperature of about 0.88 Tq. 
The mean value of these distributions is graphed as rjiq 
in Fig. [T] 

The mean nucleation or mean first-passage time, 
Txti, shows expected non-monotonic temperature 
dependence.'^ At higher temperatures, nucleation rates 
increase upon cooling because the barrier to nucleation 
decreases in size. In contrast, at lower temperatures, 
the process of crystallization is slowed by the onset of 
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FIG. 7. Time scales of tlie supercooled liquid mW water at 
1 bar. Computed structural relaxation times, riiq, are shown 
with blue points. Computed crystallization times, Txti, are 
shown with red points. Statistical uncertainties are smaller 
than the symbols. Dashed lines are drawn as guides to the 
eye. The grey region is where the liquid is unstable. 



glassy dynamics. At conditions where the amorphous 
phase is unstable, Txti and T^q are within two orders 
of magnitude of each other. These are average times 
for crystal nucleation and liquid structural relaxation, 
respectively. The distributions of these times are broad, 
with widths at least as large as the mean values. 
Therefore, the distributions of possible times for these 
respective processes will overlap, which is why a liquid 
state is no longer physically realizable in this region of 
the phase diagram. 

The boundary to unstable amorphous behavior is often 
referred to as the "homogeneous nucleation line." This 
terminology is possibly confusing because no significant 
barrier to nucleation exists in no-man's land. Indeed, 
studying the same model with straightforward molecular 
dynamics in the region of no-man's land, at T/Tq = 0.72, 
Moore and Molinero conclude that the critical nucleus is 
less than 10 molecules.l^ Coarsening times for relaxing 
defects in the crystal are necessarily longer than Txti ; and 
for mW water these times seem to be at least two orders 
of magnitude larger .13^ 

The most important point to take from Fig. [7] is that 
time scales for forming crystals are many orders of mag- 
nitude larger than those to equilibrate the liquid at stan- 
dard conditions. This fact explains why straightforward 
molecular dynamics simulation has thus far proved to be 
an unreliable probe of phase transitions in supercooled 
water and related materials. In the case of amorphous 
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phase behavior, close to or within no-man's land, we see 
from Fig. [7] that mW water requires time scales to equili- 
brate that are 3 orders of magnitude larger than those of 
the normal temperature liquid. Thus, while Moore and 
Molinero's study of freezing at such conditions'^^ is suf- 
ficiently long to illustrate likely trajectories leading to a 
crystal, it is too short to provide quantitative information 
on the underlying probability distributions that dictate 
the instability of the liquid phase. It is also not possible 
from that study to determine if other trajectories exist 
that might lead to a second metastable liquid phase. 

For real water (or more elaborate atomistic models) 
equilibration times must account for reorganization that 
overcomes donor-acceptor asymmetry of hydrogen bond- 
ing, a feature that is absent in mW water. Barriers to 
orientational reorganization will be comparable to those 
of translational reorganization. Thus, near the no-man's 
land boundary one expects equilibration times of, say, 
ST2 water, to be several orders of magnitude longer than 
those of mW water. Moore and MolinercPSl estimate it 
to be seven orders of magnitude longer. This is an issue 
that can be examined in future studies using methods 
of rare-event sampling.'^ For now, however, this paper 
has demonstrated that time-scale issues do not prohibit 
the systematic study of reversible phase behavior of wa- 
ter and related systems using the methods of free energy 
sampling, ^^ ®-' and such study draws a picture contrary 
to a widely popularized n otion of a secon d critical point 
at supercooled conditions .lilSlIlIIlIllIIMTl 
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FIG. 8. The restricted contracted free energy function, 
-F(p;Qjr'"'), for the mST2 model with iV = 216 at the tem- 
perature and pressure 235K and 2.2 kbar, respectively, and 
A/i = 0. The orientational order parameter Qg is restricted 
to Qa ^ Q™"" = 0.13. The left panel does the integration 
specified in Eq. (12) with the data presented in Section IV B, 
specifically Fig. 6. The right panel does the integration with 
data obtained with further sampling. Error bars indicate one 
standard deviation for the sampled data. The black lines are 
guides to the eye drawn through the mean values obtained 
from each sampling bin. 



tions about our results for the mST2 model. Spurred by 
this interest, this Supplement provides additional infor- 
mation. 

The results we present now focus on a contracted free 
energy function that is restricted to a range of Qg'Values, 
i.e.. 



.-F{p:.Qr^)/k^T ^ 



Jo 



dQe e 



-F{p,Qe)/kBT 



(12) 



Simulation estimates for this function when Qq is re- 
stricted to liquid-basin values are shown in Fig. [8] The 
graph on the left, obtained from the same data used to 
construct Fig. 6, shows statistical uncertainties of the 
order of ksT. The graph on the right, obtained from fur- 
ther sampling, shows uncertainties that are a fraction of 
that size. The graph on the left cannot exclude the pos- 
sibility of a subtle shoulder on the low density side of the 
density distribution. Such a shoulder could conceivably 
lead to a bi-stability at lower pressure, with a barrier be- 
tween the two basins no larger than k^T. This behavior 
is what has been reported in two simulation studies of 
the ST2 model. ■^^ ^^" For this range of densities and Qe, 
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FIG. 9. To within additive constants, the left graph shows 
the restricted contracted free energy function, F{p; Qg"*"), for 
A?^ = 216 at several indicated choices of Q™^''. These functions 
are obtained from the integration in Eq. (12) with the unre- 
stricted free energy function shown at the bottom right. The 
thermodynamic state point for this unrestricted free energy 
is (235K, 2.2 kbar) and A^ = 0.27kBT. Error bars indicate 
one standard deviation. Also shown, upper right, is the mean 
value of Qfj as a function of the maximum order-parameter 
value. Dashed vertical lines in the graphs on the right are at 
the restricting maximum values for Qg in the graphs shown 
on the left. 



however, the graph on the right excludes the possibihty 
of bi-stability with barrier heights and basin locations 
reported in Refs. 1551 and [5^ 

Dynamics of supercooled water can be glassy, and in 
that circumstance, only a limited range of Qe values will 
be explored by finite-time trajectories without impor- 
tance sampling. As such, it is interesting to consider the 
effect of restricting the range of Qe- These effects may 
explain behaviors found by others studying this model. 
Figure [9] shows its effect on the contracted free energy 
of mST2 water. The particular condition considered is 
where the free energies of the liquid and crystal basins are 
equal. Restricting Qq to values less than 0.4, however, 
does not allow the system to reach the crystal basin, and 
the average of Qe is an order-of-magnitude smaller. At 
the same time, the free energy as a function of density 
exhibits an inflection or slight minimum. This feature 
could be confused with a second liquid basin, but in fact, 
it is due to the barrier separating liquid from crystal. 
The rapid increase in (Qg) as Q^^^ increases from 0.4 to 
0.5 is indicative of the barrier that separates the amor- 
phous and crystal basins. It is only by sampling the full 
range of Qg, that the free energy shows the phases to be 
in coexistence and the mean value of Qe to be indicative 
of a crystal in a finite simulation system. 
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